Quantum Phase Transitions in Superconducting Arrays 
with General Capacitance Matrices 
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We investigate quantum phase transitions in two-dimensional superconducting arrays with gen- 
eral capacitance matrices and discrete charge states. We use the perturbation theory together 
with the simulated annealing method to obtain the zero-temperature phase diagrams, which dis- 
play various lobe-like structures of insulating solid phases, and examine the possibility of supersolid 
phase. At nonzero temperatures, an effective classical Hamiltonian is obtained through the use 
of the variational method in the path-integral formalism, and the corresponding phase diagram is 
found approximately. The insulating lobes of the solid phases are shown to exist at sufficiently low 
temperatures, and results of Monte Carlo simulations are also presented. 
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I. INTRODUCTION 



Two-dimensional (2D) superconducting arrays, weakly 
coupled by Josephson junctions, have attracted much at- 
tention [y . In the classical array, where the charging en- 
ergy is neglected, the relevant variable is the phase of the 
superconducting order parameter at each grain, and the 
logarithmic interaction between vortices is well known to 
lead to the Berezinskii-Kosterlitz-Thouless (BKT) tran- 
sition Recent development of the fabrication tech- 
niques, on the other hand, allows to make regular ar- 
rays composed of very small superconducting grains. In 
such arrays the charging energy cannot be neglected 
any longer [^, and studies considering both the self- 
capacitance Co and the junction capacitance Ci are re- 
quired. In the presence of Ohmic dissipation, the charges, 
which are conjugate to the phase variables, take con- 
tinuous values and the corresponding phase diagrams 
have been obtained through the use of a variational 
method This has revealed the existence of a low- 
temperature reentrant transition, which appears consis- 
tent with the recent quantum Monte Carlo study 
Without dissipation, in contrast, the charges change in 
discrete quanta 2e with — e being the electron charge. In 
this case the phase diagrams have been investigated via 
the coarse-graining approximation, which is mean-field- 
like in nature In general, the mean-field approxi- 
mation does not give reliable results in two dimensions, 
where fluctuations are too strong to be neglected. It is 
thus desirable to study the 2D superconducting arrays, 
with general capacitance matrices and discrete charge 
states, beyond the mean-field approximation. 

The quantum phase transitions in the superconducting 
arrays have drawn much interest also in relation to the 
Bose-Hubbard model, which describes strongly interact- 
ing bosons under the competition between the kinetic- 
energy and the potential-energy eS'ects. Here the ki- 



netic energy and the potential energy correspond to the 
Josephson energy and the charging energy in the super- 
conducting array, respectively. In the presence of the 
on-site potential, the system at zero temperature ex- 
hibits a transition between the Mott insulating phase 
and the superfluid phase, displaying lobe-like structures 
in the resulting phase diagram ||^,^. In addition, nearest- 
neighbor interactions and next nearest-neighbor interac- 
tions produce various phases such as the checker-board- 
type solid, the striped solid, and the supersolid. In par- 
ticular, the supersolid phase, which is characterized by 
superfiuidity in the presence of underlying commensura- 
bility, has been studied extensively |]9|,po|. 

This paper investigates the phase transitions in 2D 
superconducting arrays with emphasis on the effects of 
charge frustration and on the competition between the 
self- and junction capacitances. At zero temperature, we 
use the perturbation expansion and the simulated anneal- 
ing method to obtain various insulating lobes. Here ra- 
tional charge frustration introduces commensurability ef- 
fects to the array, leading to the solid phases with various 
charge densities. The detailed phase boundaries depend 
crucially on the ratio of the junction to self-capacitances, 
Ci/Cq: For example, the insulating lobes become nar- 
rower and their central positions approach the charge 
frustration values, as Ci/Cq is increased. Further, as the 
system size is increased, more and more insulating lobes 
are expected to be observed in the phase diagram, sug- 
gesting the interesting possibility of the lobe at every ra- 
tional frustration in the thermodynamic limit. To check 
the robustness of the lobes against finite temperatures, 
we apply the Monte Carlo (MC) method to the system 
without the Josephson term and find the persistence of 
the insulating lobes with low-order rational values of the 
charge density. We also use a variational method to ob- 
tain an effective classical Hamiltonian describing the sys- 
tem at finite temperatures. The corresponding phase di- 



agrams are found approximately in a self-consistent man- 
ner, which again display lobes with the half-filled charge 
density. In the absence of charge frustration, the effective 
classical Hamiltonian has only real terms, which allows us 
to use the MC method: We obtain the critical strengths 
of the Josephson coupling and compare them with ones 
obtained from the self-consistent approximation. In the 
self-charging limit, the possibility of a reentrant transi- 
tion is also pointed out. 

This paper is organized as follows: Section || intro- 
duces the quantum mechanical Hamiltonian for the su- 
perconducting arrays with general capacitance matrices 
and discrete charge states. Such a superconducting ar- 
ray is described by quantum mechanical conjugate vari- 
ables, charges and phases, and charge frustration is in- 
duced by applying an external voltage between the ar- 
ray and the substrate. In Sec. [II, we study the zero- 



the nature as well as the existence of the supersolid 
phase is discussed. Section IV is devoted to the inves- 



temperature phase transitions in the system via the per- 
turbation expansion together with the simulated anneal- 
ing method. The phase transitions at vanishing charge 
density is investigated in Sec. Ill A, where the results 
in the nearest-neighbor charging limit are al so com pared 
with those of the previous works. In Sec. IIIB| , phase 
transitions at finite charge densities are studied, and 



tigation of the finite-temperature phase transitions by 
means of a variational method and MC simulations. The 
system without the Josephson coupling is described by 
only charge variables, a nd w e present the correspond- 
ing MC results in Sec. [VA, where the charge densi- 
ties are obtained as functions of the charge frustration 
at various temperatures. The resulting step structures 
at zero temperature are compared with those obtained 
in Sec. 



[II B 



Whereas the phase diagrams in the pres- 
ence of the Josephson coupling are obtained in a self- 
consistent manner in Sec. IV B, the comparison with MC 



results reveals that the self-consistent approximation is 
accurate only at sufhciently high temperatures. Finally, 
a summary of the paper is given in Sec. |^. 



II. MODEL HAMILTONIAN 

We consider an L x L {N = L?) square array of Joseph- 
son junctions. The extra charge Qi on the superconduct- 
ing grain at site i can be written as 



^Ecu*.' (1) 

j 

where Co and Ci are the self- and the junction capacitances, respectively, and $i is the electric potential of the 
grain i with respect to the substrate. If we apply an external voltage between the substrate and the ground, 
Eq. (|l|) is changed to Qi — ^ - Cij^j+Qx, with the "gauge charge" Qx = Cq^x- The Hamiltonian describing such a 
superconducting array consists of two parts: 

H = Ho + V 

= \ - Q-)C~j'{Q, - Q.) ~ Ej E cos((^, - (2) 

i,J (ij) 



where 0i represents the phase of the superconducting or- 
der parameter at site i. Note that the charge Qi and the 
phase (pi are quantum mechanical conjugate variables sat- 
isfying the commutation relation [(f>i,Qi] = 2ei (e > 0). 
In the charging energy term Hq, the charges at sites i 
and j interact via the inverse of the capacitance ma- 
trix C~j^ whereas the summation in the Josephson en- 
ergy term V is to be done over all the nearest neigh- 
boring pairs. When the lattice constant of the array ap- 
proaches zero, we can write Eq. (|l|) in the continuum form 
g(r) = Co$(r) - CiV^^{r) and obtain $(r) due to a 
single positive charge at the origin ||ll[: <I>(r) oc Ko{r/A) 
with the modified Bessel function Kq^ where the screen- 
ing length A (= y/Ci/Co) measures the interaction range 
between charges. The above Hamiltonian in Eq. (0) is 



closely related to the Bose-Hubbard Hamiltonian and to 
the spin- 1/2 XX Z Hamihonian 0. The char gmg en- 
ergy term in Eq. drives the system to arrange its 
discrete charges such that the average charge per site 
is as close to Qx as possible, and the Josephson energy 
corresponds to the hopping energy in the Bose-Hubbard 
model . Consequently, if the charging energy is suffi- 
ciently larger than the Josephson energy, hopping is sup- 
pressed and all the charges are arranged to form a peri- 
odic lattice, leading to the insulating phase at zero tem- 
perature. In particular, when Ej — 0, the energy eigen- 
state can be chosen also as a simultaneous eigenstate of 
the charge operator since [Qi,H] = 0. The uncertainty 
relation between the charge and the phase then asserts 
that phase coherence and consequently, superconductiv- 



2 



ity cannot be attained in this limit. In the opposite limit 
of the strong Josephson coupling, all the charges are in 
the extended state owing to the dominant hopping term, 



and the superconducting phase is favored. For conve- 
nience, we write Qi = Qi/2e, Qx = Qx/2e, Eq = e^/2Co, 
and Cij = Cij/Co, and obtain Eq. (0) in the form 



H = Ha + V = AEoJ2il^ 



^,3 



Ix) - cos((/)j 



(3) 



where the inverse of the dimensionless capacitance ma- 
trix is given by 



1 



gik-(xi-Xj) 



N ^1 + 4Ci/Co - 2(Ci/Co)(cos + cos ky) ' 



(4) 



III. ZERO-TEMPERATURE PHASE DIAGRAMS 

In general the state which minimizes the energy deter- 
mines the zero-temperature phase of the system. Sup- 
pose that we have obtained the energy levels of the sys- 
tem which correspond to the ground state, the first ex- 
cited state, etc. If we introduce a small change in the 
control parameters, the energy levels shift by a small 
amount, and it is possible that the ground state and the 
first excited state cross each other as we pass through a 
set of values of the control parameters; this "level cross- 
ing" leads to a zero-temperature phase transition in the 
parameter space [|3j. It should be noted here that con- 
sideration of only the two low-lying states, the ground 
state and the first excited one, may not be sufficient. 
Figure |i] shows schematic dependence of the energy levels 
En on the control parameter K in two possible cases. In 
Fig. H (a), the ground state and the first excited state 
cross at the critical value Ki of the control parame- 
ter, which corresponds to the true transition point. In 
Fig. |l| (b), on the other hand, the true transition point 
is given by K2, where the ground state and the second 
excited one cross each other. In this case consideration 
of only the two low-lying states would lead to Ki , which 
is not the true transition point. To find the true transi- 
tion point, we thus need to consider all the energy levels, 
which is not possible in practice. Despite this, consid- 
ering the level crossing of the ground state and the first 
excited state is very useful and convenient to obtain the 
upper limit of the true transition point. 



In terms of charges the system can be either in the 
incommensurate conducting phase or in the commen- 
surate insulating phase; in the latter charges are ex- 
pected to form a superlattice which covers periodically 
the whole system. If the commensurate state has en- 
ergy lower/higher than that of the incommensurate state 
for given values of parameters, the system should be in 
the insulating/conducting phase. (Here Ej/Eq, q^, and 
Ci/Cq constitute the control parameters.) Accordingly, 
for given total number of charges, it is needed to find the 
charge configurations of the commensurate state and of 
the low-lying incommensurate states. The symmetry of 
the Hamiltonian under the transformations 



Qx + 1, 
-Qx, 



(5) 



allows us to consider only in the range < < 1/2. 
In this restricted range < g^; < 1/2, the commensurate 
states of the system may be classified in terms of the av- 
erage number of charges, (qi) = J2i Qi/^- For (qi) = p/q 
with relatively prime integers p and g, the charges are 
expected to form a q x q superlattice, on the unit cell of 
which pq charges, each with the unit charge, are located. 



A. Phase transitions for (qi) — 

We first investigate the phase transitions in the system 
with (qi) ~ 0. The energy of the (g^) = commensurate 
state is computed via the second-order perturbation ex- 
pansion with the Josephson term as a perturbation in 
Eq. (0), and is written as 



i?(,)(fe) = o) = i?g(o) + £;W(o) 



(6) 



Without the perturbation, the Hamiltonian of the system 
contains only charge operators and yields the eigenfunc- 
tion 



(0|q) ^ 



&jv|gi,g2,g3, • • • ,gw; = 



N 



exp 



(7) 
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which is simply the plane-wave state of the usual free 
particle system, with the interpretation of (pi and qi 
as the position and the momentum, respectively. For 
the commensurate state with (g^) = 0, all the charge 
eigenvalues g^'s are zero and it is easy to compute the 
zeroth-order contribution E^P} (0) and the first-order con- 



tribution i?g(0): EfJ^{Q) = ^EmlT.^.,C;^ = ^E^qlN 



and E'^^liQ) = (q = 0|y|q = 0) = (1/2^)^ / 



^-^ cos(0i — (pj) = 0. The second-order contribution is 
given by 



(q'l^|q = 0)| = 



q'#0 



£;(O)(q = 0)-£;(O)(q')^ 



(8) 



where E'^^^ (q) is the zeroth-order energy for the charge 
configuration q, and 



(q'|F|q = 0) = (q'l - i?j ^ cos(0, - </),)|q = 0) 
Ej 



2(27r) 



N 



E 



(9) 



Inserting Eq. (||) into Eq. (H), we obtain the energy of 
the {qi) — Q commensurate state to the 2nd order in Ej: 



E^^^iO) ^ AEoNql 



EjN 



8E, 



[CqJ - C^o^) 



(10) 



We next consider the incommensurate state with 
(qi) = 0. Among the various excitations present in the 
system with the constraint (qi) = 0, the lowest excitation 
is expected to be point-like and we consider two types of 
excitations: One is the creation of a single charge (SC) 
at one site and the other is the creation of a charge 
dipole (CD) at one bond [qi — 1, qj = —1 with (i,j) 
being the nearest-neighboring pair]. Both excitations 
do not change the value of (qi) in the thermodynamic 
limit {N — + oo). In the self-charging limit (Ci — 0), 
the interaction matrix C~j^ takes a diagonal form and 
the first excited state should be of the SC type, since 
the CD type excitation requires the creation energy of 



two charges. On the other hand, in the nearest-neighbor 
(NN) charging limit (Cq ~ 0), only the excitations which 
does not change the total number of charges are allowed 
and the first excited state is expected to be of the CD 
type. Accordingly, as the value of Ci/Cq is increased, 
the first excited state should change from the SC type to 
the CD type at a critical value of Ci/Cq. In particular, 
the system in the NN charging limit displays the charge- 
anticharge unbinding transition [|ll|,^ 15 1, the BKT na- 
ture of which is manifested by the square-root cusp in 
the resistance This charge BKT transition is sup- 

pressed by the presence of the self-capacitance, making 
the screening length finite, and it seems likely that the 
change of the first excited state according to the value of 
Ci/Cq in turn leads to the change in the nature of the 
transition. 

Between the above two types of excitations, the lowest 
one is determined from the comparison of the zeroth- 
order energies of both types: If AE given by 



A£; = £:(°)(sc) (CD) 

= iEo (CoV - 2q, + Nql^ - AEq {c-; - - C-; + C-,' + Nql) 



AEq 2C 



C, 



00 



(11) 



is negative/positive, the creation of a single charge/charge dipole is the lowest excitation. Figure^ shows the boundary 
for AS = in the limit of ^ cxd. It is found that for Ci/Cq < 14.116 the lowest excitation is of the SC type 
regardless of the value of qx, and that the region of the CD type excitation (the region below the curve in Fig. ^) is 
quite small, allowed only for q.^ < 9.4 x 10""*. Therefore we below focus on the SC type excitation for Co 7^ and 
qx > 0, and then consider the CD type excitation only in the NN charging limit (Co = 0) with q^ = 0. 
We write the energy of the first excited state of the SC type as 



i?(,,)((g.> = 0) E^^^(0) + <\ (0) + <\(0) + 0(i?5), 



(12) 
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and the eigenstate as 



|q(')) = kf\ (J)) with (13) 

' = otherwise, 

which is expected to be a superconducting state since the extra charge at site I can hop to any site without energy cost. 
Inserting Eq. (|l|) to Hq in Eq. (|), we obtain the zeroth-order contribution: i^l^j) (0) — 4:Eo{Cqq^ — 2qx + Nq'^), where 

we have used the translational symmetry C^-^ = C^J and J2j — 1- Since the above excited state in Eq. ( p^ ) 
has an A^-fold degeneracy {I can take values 1, 2, • • • , iV), we use the degenerate perturbation theory to calculate the 
first-order contribution e'^^1^{Q). We thus need to diagonalize the matrix V whose components are given by 

-Ejjl for (/, m) nearest neighboring pair, , > 

otherwise. ^ ' 

For that purpose, it is convenient to use the Fourier transform: 

^™ = ^E^''''^'"""'"^^k, (15) 

k 

which gives straightforwardly the eigenvalues Vk ~ ~Ej(co8 + cos ky). The first-order contribution is then given 
by the lowest eigenvalue: ^^(^^^^^(O) = Vk=o = —'2Ej. For this lowest-lying excited state with k = 0, it is also 
straightforward to obtain the second-order contribution, which leads to the energy of the incommensurate state 

(0) « 4Eo{Nql ~ 2q, + CoV) - 2E.J + — ^ 



E^j 2 1 

— -r^^ ^-J ^-j- + -r^—^ T^r—^ r^r—^ ^-j" 

8^^0 \C^Q ~ ^i+y,0 ~ ^00 +^"10 C'lO ~ ^2i!ifi~ ^00 +^"10 



(16) 



Here denotes the summation over i and its four nearest-neighbors j with the restriction j 7^ 0. The zero- 

temperature transition occurs when the energy of the commensurate state given by Eq. ^ and that of the incom- 
mensurate state given by Eq. (^6|) becomes equal. This reveals the zero-temperature phase transition in the system 
with {qi) — 0, as qx is varied. Namely, the system becomes superconducting as qx is increased beyond the critical 
value g^, which is given by 



'^^ 2 °° 2E, ) 256 Uo M C„V - C^.' ^^^0 ^io' Cr^' + C^J 



I ( EjV ( 2 1 



64 \ i?o / Y '^i+y.O ^00^ + ^2A,0 ^ ^00^ + ^^10^ . 



(17) 



The resulting phase boundary between the insulating phase and the superconducting phase is shown in Fig. ^, obtained 
for a 24 X 24 array. The junction capacitance here tends to enhance superconductivity, which is in accord with the 
result in Ref. Q. It is also obvious that as the system approaches the NN charging limit, the insulating region shrinks 
to zero; this apparently cures the well-known artifact of the coarse-graining approximation, where the system remains 
in the insulating phase for arbitrarily large values of Ej in the NN charging limit Q . 

To check the validity of the second-order perturbation expansion, we have computed the third-order contribution 
in the self-charging limit, and obtain the critical value q^' 

'-<'--°'^5-Ki)-Ti(iy-T^(i)'- <'^' 
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The three curves in Fig. ^ represent the phase boundaries 
obtained from Eq. ([1^) up to the first-, the second-, and 
the third-order contributions, respectively. It is shown 
that the third-order contribution does not change sig- 
nificantly the boundary obtained from the second-order 
calculation, inducing 8.3% of the difference in the critical 
value of Ej/Eq at Qx =0. 

We next consider the NN charging limit in the absence 
of charge frustration {q^ = 0). In this case, Eq. (10) gives 
the energy of the commensurate state in the form 



E, 



EjN 



8Ei (Cqo^ - C-i^q) ' 



where Ei = e^/2Ci and 



4 — 2(cos kjc + cos ky) 



(19) 



(20) 



In the thermodynamic limit, we replace the discrete sum 
by the integral and obtain 



C, 



00 



1 



10 



(2^)2 
1 



dk 



1 — cos fcr 



2tt 



dk,, 



4 — 2(cos kx + cos ky) 
1 — cos fcj, 1 
3 — cos kr 4 ' 



(21) 



leading to ~ —EjN/2Ei. The first-excited in- 

commensurate state in the NN charging limit is of 
the CD type (see Fig. ^ with the zeroth-order energy 

E^°\=8Eii^oo'-^iQ') =2-E^i and the null first-order en- 



ergy E 



{ic) 



0. It is complicated but straightforward to 



(2) 

calculate E^^^^ from the second-order degenerate pertur- 
bation theory, which leads to a 4N x 4N matrix. The 
Fourier transform then yields TV 4 x 4 Hcrmitian ma- 
trices, each of which can be diagonalized. The lowest 
eigenvalue then gives the energy of the CD type incom- 
mensurate state, up to the second order: 



E, 



2E, 



8Ei 



(47V 62.000) , 



(22) 



where the number 62.000 has been obtained numerically 
for a sufficiently large array (1024 x 1024); at this size 
finite-size effects have been confirmed to be negligible. 
The critical value of Ej/Ei, below which the array is 
in the insulating phase, is obtained from the condition 



E{ic) - E(^c) = 0: 



Ej , 

-f I « 0.508. 
El 



(23) 



Note that this value is somewhat larger than the criti- 
cal value {Ej/Ei)c « 0.23 obtained in Refs. ||,|ll|. One 
possible explanation of the above discrepancy is the fail- 
ure of the second-order perturbation expansion. Namely, 
consideration of higher-order terms may eliminate the 
discrepancy. The other possibility is that the CD type 
excitation does not give the true transition point: See 
Fig. (b) for a schematic picture, where Ki may corre- 
spond to the obtained transition point {Ej/Ei)c « 0.508 
and K2 the true transition point. In experiment, the 
even larger critical value {Ej/Ei)c « 0.6 has been ob- 
tained [0 , via the curve-fit to the square- root cusp form 



of the resistance [Q. Recent quantum Monte Carlo 
simulations, on the other hand, have yielded the result 
{Ej/Ei)c = 0.36 ±0.04. 



B. Phase transitions for (qi) 7^ 

In the NN charging limit, there exists an interesting 
duality between the charge and the vortex |11|, which 
concludes that in the absence of the Josephson term, the 
charges form a superlattice determined by the value of 
charge frustration q^. This is a counterpart of the vor- 
tex superlattice in the classical array without the charg- 
ing energy, formed in the presence of magnetic frustra- 
tion 17|. In the system considered here, on the other 
hand, the non-vanishing self-capacitance gives the inter- 
action range between charges A — \/Ci/Cp in units of 
the lattice constant, thus making it finite [Q. Accord- 
ingly, we expect that the formation of the q x q charge 
superlattice with g ^ A is improbable in the system with 
the rational charge frustration q^ — p/q. This implies 
that (qi) is not necessarily equal to qx in the array with 
non- vanishing Co, which reflects that a short-ranged po- 
tential in general makes the commensurate phase extend 
over a finite range of the chemical potential. (Note that 
the charge frustration in this work plays the role of the 
chemical potential.) 

We now investigate the phase transitions for (qi) = p/q 
(7^ 0), using the canonical ensemble with the total num- 
ber of charges fixed. To find the zero-temperature phase 
diagrams we need to know the configurations of the com- 
mensurate states and of the lowest-lying incommensurate 
states, which is a highly nontrivial problem in the 2D 
system considered here, especially for large q. For the 
ID Josephson junction array with (q^) — 1/q, the com- 
mensurate state, the first excited state, and the resulting 
zero-temperature phase boundaries have been obtained 
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via the perturbation expansion |18|. We write the en- where the relations J^i^i — ^{Qi) = ^p/q a-nd 

ergy of the commensurate state with (qj) = p/q 7^ in Cjj^ 1 give Ej^^ip/q) = ^EoY^QiCij^l] - 

the form 8Eoq^Np/q + AE^qlN . In general, when {q,) ^ 0, 

E{c){p/q) = E\^J^{p/q) + E\^J^{p/q) + E)^^!^{p/q) + 0{E^j), the commensurate ground states are degenerate, which 



(24) 



makes it necessary to consider the matrix element: 



n*(e 



(m) 

9fe 



(m) 



r / (0 (m) ,\ . / (0 
, r / (0 l™) I 1^ X / (0 (™) 



+ 1 



(25) 



with |q^'^) being the Ith. degenerate commensurate state. Note that the configuration of the superlattice unit cell in 
a degenerate commensurate state is in general different from that in a different (degenerate) state. Therefore, for 
I ^ m, q('^ differs from q^™) at the total of 0{N) sites for finite q, leading to Vim = 0. Further, it is obvious that 
Vu = 0. It is thus concluded that all the components of the matrix V vanish, which results in E^^^{p/q) = 0. The 
second-order contribution is given by 



(26) 



where i?*^°)(q) and E^^'^lq') is the zeroth-order energy for given charge configuration q and q', respectively, and the 
summation is over the states q' making the denominator nonzero. The superscript I on q has been dropped 

(2) 

since all the degenerate commensurate states give the same contribution to E^^J^ip/q). Namely, the degeneracy is 

not removed in a finite order since q^'^ and q(™) cannot be connected by a finite number of charge hoppings. With 
Eq. (p5|), Eq. (p6|) reads 



E^/A 



i?(0)(q)-i?(")(q(,.,))' 



(27) 



where <\{ij) = {qi,q2, ■ ■ ■ , qi-i, (Zi+l, qi+i ■ ■ • , q_j-i, • • ■ 7 Qn) is obtained from q by a charge hopping between 

sites i and j. It is easy to compute the denominator in Eq. (27): 

£;(0)(q) _ ij(0)(q(,^^.)) = AEoJ^ [(% - <l--)Cki\li - - (^(^j),^ - q.)Ci:i\q^,,,^,i - q.) 



k,l 



8E0 



•-"10 



(28) 



which leads to the expression for the energy of the commensurate state with (qi) = p/q: 



E/ 
'32Eo 



Y.' 



^-1 



qi 



^00 ~ ^10 



(29) 



Here qi is the charge at site i in the commensurate state with (qi) = p/q. Since the (commensurate) ground state 
configurations are not known for general p/q-, qi in this work is obtained via the simulated annealing method. 

To compute the energies of the incommensurate states with {qi) ^ p/q 7^ 0, we classify excitations according to the 
number q^ of extra charges (with the sign taken into account), and let E(^icj{p/q, q^) denote the energy of this excited 
state which has the total number of charges Np/q + q^. To the second order in Ej, it is straightforward to obtain the 
energy of the gg-charge excited state: 



7 



E 



32Eo 



1 



(ij) Y.I ipji" - ^) 9^'^''' ~ (Cqo^ - C'lo^) 



(30) 



where is the charge at site i in the gg-charge ex- 
cited state, and the first-order contribution E^f^\^{p/q^qe) 
has been observed to vanish in numerical calculations, 
regardless of the values of p/gi^ 0) and q^ considered 
in this work. Further, in Eq. (pQ), we have assumed that 
degeneracy is not removed in the second order although 
there are some exceptions, e.g., for pjq— 1/2 and qe = —1 
(see below). To find the phase boundary, we need to 
know the charge configuration which corresponds to the 
first excited state for given values of Ci/Cq and qx- For 
this purpose, we use the following numerical procedure: 

1. Distribute Np/q + q^ charges randomly. 

2. Use the simulated annealing method to find 
the charge configuration ql'^''\ which minimizes 

the zeroth-order energy A^EoJ^ijli'^'^^ij^lj'^'^ ~ 
8Eoqx{Np/q + qe). 

3. Change the value of ge, and repeat the steps 1-3. 

4. Among the obtained charge configurations with 
their zeroth-order energies, find one corresponding 
to the lowest-lying incommensurate state. 



In the numerical simulations, we have used the con- 
ventional simulated annealing method adopting the 
Metropolis algorithm for the arrays of sizes 10 x 10 and 
12 X 12. Starting from high temperatures (T > 1.0), 
we decrease the temperature slowly with the decrements 
AT = 0.1 (at high temperatures) and 0.002 (at low tem- 
peratures), and perform 10000 MC steps per site at each 
temperature step. We have found that the charge config- 
uration corresponding to the lowest energy depends cru- 
cially on the value of Ci/Co, and thus used three different 
values: Ci/Cq = 0.1, 1.0, and 5.0. Through 20 indepen- 
dent runs, we have obtained the charge configurations 
corresponding to the lowest energy, which have also been 
checked via the entropic annealing algorithm |l^. Al- 
though we have allowed double occupancy {qi = 2) as 
well, in the lowest-lying configurations qi is always found 
to have the value either or 1. 

Once the charge configuration gj^"'* which minimizes 
the zeroth-order energy in Eq. (^0|) is found, the phase 
boundary is determined from the condition E(^^-j(j}/q) = 
E{ic){p/q, 9e)- This leads to the critical value q^ separat- 
ing the superconducting phase from the insulating phase: 



2qe 



+ 



{Ej/E^f 
256(7e 

-E' 
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where {q^} again describes the charge configuration of the 
commensurate state. When q^ < 0, the system, display- 
ing superconductivity for q^ < q^, becomes insulating as 
q^ is increased beyond gj. This behavior according to 
whether qx is larger or smaller than q'^ is reversed when 
qe > 0. 

Table | shows the values of qe in the lowest excitations, 
breaking the commensurability of the states with {qi) = 
1/2, 1/3, 1/4, and 1/5, in the arrays of sizes L = 10 and 
12. Except for {q/) = 1/2, there are two values of qe in 



each case, one for the positive-charge excitation {qe > 0) 
and the other for the negative-charge excitation {qe < 0). 
Since the value of the charge frustration has been re- 
stricted in the range < < 1/2, we need to consider 
only the negative-charge excitations for (qi) = 1/2. It 
displays that as Ci is increased, the absolute value of qe 
decreases, manifesting that the excitations become more 
point-like. As an example, we present the charge config- 
urations for (qi) ~ 1/2 in Fig. |[ where (a), (b), and (c) 
correspond to the commensurate state, the excited state 
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with = —1, and the excited state with = —12, re- 
spectively, in a 12 X 12 array. It is of interest to note that 
the configuration (c), which has 60 charges, is precisely 
the same as the ground state configuration of vortices in 
the presence of the magnetic frustration 5/12 When 
(7e = — 1 is the lowest excitation, there exists A^/2-fold 
degeneracy since the vacant site, represented by a small 



in place of Eq. (^?0|). Here k denotes the position of the 
vacant site. 

Figure ^ shows the resulting phase diagrams for 
Ci/Cq= (a) 0.1, (b) 1.0, and (c) 5.0. The various lobe- 
like structures, in which the array is in insulating phase, 
are obtained. Since the screening length A = Ci/Co, 
beyond which the magnitude of the interaction between 
charges decreases exponentially, is much smaller than the 
system size, the finite-size effects on the phase boundaries 
are presumably negligible. We have checked the phase 
boundary for {qi) — 1/2, and indeed found that it hardly 
changes with the system size considered here (see below) . 
Accordingly, the phase boundaries in Fig. ^ are expected 
to be qualitatively correct even in the thermodynamic 
limit. As we increase the array size, more insulating- 
lobes should be observed, since the array with Ej = 
should be insulating, regardless of the value of qx- At 
this stage, however, it is still unclear whether the insu- 
lating lobe exists at every rational value of q^. In the NN 
charging limit, the charging-energy term takes the form 
Hq = 4:Ei J2t,jiQt - 1x)C~^{qj - Qx) with C^^ given by 
Eq. (pO|). Here the divergence of Cqq leads to (qt) — qx, 
which is simply the counterpart of the vortex-neutrality 
condition in the classical array. Indeed Fig. || shows that 
the position of the insulating region approaches (qi) = qx, 
as Ci/Co is increased. 

To study the effects of the junction capacitance versus 
the self-capacitance in detail, we now concentrate on the 
(qi) = 1/2 insulating phase, whose commensurate charge 



square in (b), can be located on any of the total N/2 
sites. In this case, the degenerate states are coupled in 
0{Ej), and we should use the second-order degenerate 
perturbation theory to calculate the energy. It is again 
straightforward to diagonalize the second-order matrix 
by means of the Fourier transform, which yields 



(32) 



configuration displays the 2x2 superlattice structure. 
[See Fig. || (a).] Similarly to the vortex superlattice in 
the fully- frustrated classical array the commensu- 
rate state has the two-fold degeneracy. If we remove 
a single positive charge from the system, commensura- 
bility is destroyed near the resulting point defect. The 
charge configuration of the corresponding lowest-lying in- 
commensurate state has been shown in Fig. |] (b), which 
describes the first excited state for Ci/Cq = 1.0 or 5.0. 
When Ci/Cq = 0.1, in contrast, the first excited state has 
the configuration shown in Fig. ^ (c), with the number of 
extra charge qe — —12. This observation that the nature 
of the lowest excited state changes with Ci/Cq has an im- 
portant implication with regard to the supersolid phase, 
where a commensurate charge-density wave is believed 
to coexist with superfluidity ^,|l^: The charge configu- 
ration in Fig. H (b) regarded as a snapshot of the charge 
configuration in the supersolid phase demonstrates that 
the underlying global commensurability still holds de- 
spite of the local point defect. Furthermore, since the 
point defect can hop to any site without energy cost, it 
is delocalized in the energy eigenstate, leading to super- 
fluidity. Accordingly, we conclude that the existence of 
the supersolid phase depends on the value of Ci/Cq. In 
particular, when Ci/Cq is sufficiently small, the super- 
solid phase does not exist. It is of interest to compare 
this argument pointing out the importance of the vacant 
defect in supersolidity with that of Ref. , emphasizing 
the role of the double occupancy. In the numerical inves- 
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tigations to find the lowest-energy charge configurations, 
we have found that all the charges have the value either 
or 1 at zero temperature, and that those configurations 
with Qi = —1 or 2 at some sites have much higher ener- 
gies. The (qi) — 1/2 phase boundaries in Fig. (b) and 
(c) may be considered to divide the insulating commen- 
surate phase and the supersolid phase. As we go further 
into the supersolid phase starting from the phase bound- 
ary, the more point defects are generated; the commen- 
surability should be completely destroyed when we cross 
another phase boundary separating the supersolid phase 
from the superfluid phase, although within the formalism 
adopted here, we are unable to find this new boundary. 
To obtain the phase diagram for (qi) = 1/2, we have con- 
sidered L X L arrays up to L = 24, and found that the 
finite-size effects are negligible for Ci/Cq < 10, which 
corresponds to the interaction range A < \/lO <C L, and 
that the qualitative features of the phase diagram do not 
change even when Ci/Cq ^ 10. 

We have further found that the difference between the 
phase boundary obtained from the excitation q^ ~ —1 
and that from q^ — —12 is insignificant in the 12 x 12 
array with Ci/Cq — 0.1. Therefore, for simplicity, we 
consider only the excitation Oe = — 1 and compute the 
phase boundary from Eqs. (|2^) and (^2|). The resulting 
phase diagram obtained for the array of size L = 24 is 
shown in Fig. |^. The insulating region is shown to ex- 
pand with Ci/Cq when Ci/Cq ^0.1 and then to shrink 
as Ci/Cq is increased further. This is to be compared 
with the case {qi) = (shown in Fig. ||), where supercon- 
ductivity is enhanced monotonically with Ci/Cq, even 
for small Ci/Cq. 



IV. FINITE-TEMPERATURE PHASE 
DIAGRAMS 



In this section, we investigate the phase transitions of 
the superconducting arrays at finite temperatures, where 
thermal fluctuations as well as quantum fluctuations sup- 
press the ordering of phases. While at high temperatures 
the system should be in the normal (disordered) state 
due to large thermal fluctuations, strong charging effects 
at low temperatures produce large quantum fluctuations, 
tending to destroy superconductivity of the system. We 
thus expect rich behaviors resulting from the interplay 
of thermal and quantum fluctuations, which may be con- 
trolled by the temperature and charging energy. We first 
study the system in the absence of the Josephson cou- 
pling via MC simulations, and then investigate the phase 
transitions in the presence of the Josephson coupling via 
a variational method together with MC simulations. 



A. Superconducting array without Josephson 
coupling 

In the absence of the Josephson coupfing [Ej — 0), the 
Hamiltonian is diagonal in the charge representation and 
written in the form 



(33) 



where qi is the eigenvalue of the charge operator at site i. 
We perform MC simulations on the above Hamiltonian to 
investigate how the insulating phase corresponding to the 
regions inside the lobes in Fig. ^ changes with the tem- 
perature. Using the conventional MC method with the 
Metropolis algorithm, we start from a sufficiently high 
temperature (T = 1.0), and decrease the temperature 
slowly with the decrements AT = 0.1 (at high tempera- 
tures) and 0.002 (at low temperatures). During the simu- 
lations, we allow creation and annihilation of a charge at 
a site, as well as hopping to other sites. At each temper- 
ature step, we compute the thermal average {qi)T using 
10000 MC steps per site, after 1000 MC steps of eqmli- 
bration. 

Figure || displays the results for arrays of sizes L = 10 
(left) and 12 (right), with Ci/Cq = (a) 0.1, (b) 1.0, 
and (c) 5.0. At zero temperature, the curves represent- 
ing (qi) versus q^ exhibit steps at rational values giv- 
ing the (commensurate) insulating phase. Such steps 
structures, which correspond to the lobe structures in 
Fig. ^, are common in systems displaying commensurate- 
incommensurate transitions, and signal the commensu- 
rate locking |^l|. It is shown that, as expected, the width 
of a step decreases with Ci/Cq, and the curve approaches 
the straight line (qi) = qx in the limit Ci/Cq oo, man- 
ifesting the suppression of the insulating phase. Compar- 
ing the results for the L = 10 array with those for the 
L = 12 array, we also expect more rational steps in larger 
arrays, which is apparently suggestive of the devil's stair- 
case structure pl| in the thermodynamic limit. For con- 
clusive results, however, more detailed study is needed. 
Table || presents the widths of the rational steps of 
{qi) = 0, 1/3, and 1/2 in the L = 12 array at zero tem- 
perature. One can see here that the values in Table ^ 
are in good agreement with the Ej = results in Fig, 



which supports the validity of the approach in Sec. [II. 

Figure || also shows that steps tend to disappear at 
finite temperatures, reflecting that thermal fluctuations 
suppress the charge ordering due to quantum coherence 
as well as the phase ordering. Accordingly, the insulat- 
ing commensurate phase in general turns into the nor- 
mal (disordered) phase as the temperature is increased 
from zero. Note, however, that some steps survive weak 
thermal fluctuations. Namely, whereas steps at higher- 
order rationals easily disappear, low-order steps such as 
(qi) = 0, 2, and 1/3 persist at low but nonzero temper- 
atures. These features presumably can also be observed 
in the systems with nonzero Josephson couplings. Thus 
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as we increase the temperature from zero, higher-order 
insulating lobes in Fig. y disappear, but some low-order 
lobes are expected to remain at low temperatures. 



B. Superconducting array with Josephson coupling 

In the presence of the Josephson coupling, both 
the charge and the phase variables should be treated 
quantum-mechanically, which prohibits exact analyti- 
cal treatment. Here, we use the Giachetti-Tognetti- 
Feynman-Kleinert (GTFK) variational method [2^ to 
evaluate the path integral and to obtain the effective 
classical Hamiltonian, which is convenient for studying 
the phase transitions at finite temperatures. The GTFK 
method, which remains reliable in known cases even at 
zero temperature, has been successfully applied to the 
superconducting arrays with continuous charge states in 



the presence of Ohmic dissipation [p|j23| . We begin with 
the Hamiltonian given by Eq. (|^) and write the partition 
function of the system in terms of the path integral : 

Z^^llJdMO) JvMr)e^p{~SE[Mr)]} (34) 

with the Euclidean action 

Se= / drLEir), (35) 
Jo 

where P = l/fc^T is the inverse temperature, and the 
Planck constant has been set equal to unity (h = 1). 
The Euclidean Lagrangian Le can be obtained from 
the Hamiltonian via the Wick rotation t — > —it and 
the Legendre transformation. We use the representa- 
tion Qi = —id/d(f)i obeying the commutation relation 
= i, and get 



Le = 



16K 



^ ^ (f>tCij<f>j ~ iqa; ^(f>i~ Ej ^ cos( 



(36) 



with (f)i = d4>i/dT, where the charge frustration enters through a purely imaginary term. 

In the presence of Ohmic dissipation, the charge at each site takes continuous values, and only the paths satisfying 
(t>i{P) = '/'^(O) contribute to the partition function, with the dissipative action term added [^,^- In this case, the 
charge frustration does not appear in the Euclidean action since dripi — 0, and the resulting phase diagram does 
not depend on q^. If there is no Ohmic dissipation present, in contrast, only discrete charge states are allowed, and 
all the paths satisfying (/3) = (0) -I- 2'Kni with integer n,; contribute to the path integral in Eq. (^) . Therefore 
the path integral should include the summation over the "winding numbers" {rii}, yielding interesting behaviors 
associated with the charge frustration Accordingly, the allowed charge states, i.e., continuous or discrete, are 
crucial in the resulting phase diagrams. 

From the boundary condition (j)iiP) = + 27rni, we decompose 0^ into the periodic variable 9i satisfying 

9i{(3) = 9i{0) and rn, according to 0i(r) = Oi{T) + {2'kt/ (3)ni, and obtain the Euchdean action in the form 



where 



Se — 



1,3 



(37) 



Sph 



^ . ^ . r 27rT 

dr i — r X! ^ ^ + ~B~^'^' ~ ^^^^ 

Following the GTFK method, it is straightforward to obtain the effective classical Hamiltonian 



where g is determined by 
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gf3Ej cos( 
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(39) 



with xi = y^4:gof3^EjEo/{l + 2Ci/Co) and X2 = ^/igoP^EjEo/il/S + 2Ci/Co). Here we have used identity 

^2 



V I 

n—l ^ ' ' 



Yx 



2 (xcothx — 1). 



When Eq. (39) has more than one zeros, we should choose the largest one according to the extremal principle. The 
detailed procedure is entirely similar to that in Ref. [Q, and will not be repeated here. As a result, the quantum 
fluctuations renormalize the Josephson coupling Ej to gEj. Furthermore, in contrast to the continuous charge case [Q, 
the charge frustration appears explicitly in the Hamiltonian, and is expected to play an important role in the phase 
transition of the system. 

Unfortunately, the Kronecker delta Sni,nj between winding numbers of nearest neighboring pairs in Eq. ( p8| ) makes 
further analysis very difficult. For simplicity, we thus replace it by its self-consistent average S and obtain the 
approximate Hamiltonian 



- iSH'^i = -T-rrr nidjUj + i2Trq^ + g^P^j cos(e'i - 0j), 



(40) 



where 



Tr Sr, 



Tr e-'^^e 



is the ensemble average with respect to the approximate Hamiltonian H'^^. The Poisson summation formula then 
allows us to write 



6 = 









E{m,}exp 







(41) 



with q[{x) = x(5i_i — (5i,2), which manifests the absence of 
dependence on Ej. Once 5 is determined from Eq. (41), 
the partition function of the system obtains the form 



Z — ZchZph^ 

where the charge-part 



(42) 



Zch = ^ cxp 

{q^} 



-4/3£;o^(g.-g.)C^\gj-g.) 



has been obtained from the Poisson summation formula 
and the phase-part is given by 



Zph — 



cxp 



-g5pEj cos(0, - 



Unless the ratio Ci/Cq is infinite, the interaction C^^^ 
between charges is short-ranged and does not have any 



singularity. It is then reasonable to assume that the criti- 
cality of the system is governed by the phase-part, which 
implies that the system exhibits a vortex-unbinding BKT 
transition at the critical temperature given by pq] 



g5pEj « 1/0.9. 



(43) 



We first consider the self-charging limit (Ci = 0). 
In this limit, 5 can easily be evaluated since the 
terms in Eq. (|4^) factorize. The resulting phase dia- 
grams determined from Eq. ( ^3|) at temperatures T = 
0.1,0.2,0.5,0.7, and 1.0 (in units of Eq/Ub) are shown 
in Fig. ^ The region in the right-hand side of each line 
corresponds to the superconducting phase whereas the 
left region represents the normal (disordered) phase. In 
this phase diagram, we observe that the superconduct- 
ing region expands as qx is increased at all temperatures. 
(Recall that the symmetry of the system allows to re- 
strict within the range [0,1/2].) Therefore it is concluded 
that the charge frustration in general enhances supercon- 
ductivity. The qualitative features are largely similar to 
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those of the results from the coarse-graining approach ||] 
at temperatures T > 0.2. Unhke the latter, however, 
Fig. P also shows reentrant behavior at low tempera- 
tures, for small qx- At zero temperature, unfortunately, 
we obtain the unphysical result that the system is in the 
superconducting state for 1/3 < Qx l£ 2/3, even in the 
limit of the vanishing Josephson coupling. This is pre- 
sumably an artifact of the self-consistent approximation 
(SCA) replacing of the Kronecker delta by its average, 
with correlations neglected. 

When Ci/Cq ^ 0, on the other hand, the terms in 
Eq. ( p] ) do not factorize, and the average S cannot be 
evaluated analytically, which makes it inevitable to re- 
sort to a numerical method. At extremely low temper- 
atures, nii in Eq. ( ^l]) can have the value either or 1, 
which allows us to calculate S by direct summation of 2^ 
terms for small L. At high temperatures, other values 
of mi are allowed; this makes the direct summation im- 



practical and we use MC simulations. Figure |l^ shows 
the resulting phase diagram obtained for a 4 x 4 array 
with Ci/Co = 1.0 at temperature T = 0.1, 0.2, 0.3, and 
0.4. At T — 0.1, 6 has been calculated by means of 
both the direct summation with rrii — and 1 and MC 
simulations, which give results in good agreement with 
each other. At higher temperatures, S has been obtained 
via MC simulations, yielding the enhancement of super- 
conductivity by charge frustration in a similar manner 
to that in the self-charging limit. At low temperatures 
(T = 0.1), on the other hand, there exists a small lobe 
near qx = 1/2, which indicates the expansion of the insu- 
lating region due to the commensurability effects on the 
ordering of charges. 

To check the validity of the SCA adopted here, re- 
placement of (5ni,n to its self-consistent value (5, we have 
performed MC simulations with qx — 0, at which the 
Hamiltonian in Eq. (pSh has only real terms: 



l3Hci = ^^p^^ riiCijUj + gPEj ^ cos(6'i - 6'^) (5„;,„^. 



(44) 



We have considered a 16 x 16 array with various values of Ej/Eq, and measured the helicity modulus T with 100000 
MC steps: 



T = — 

N 



x% cos( 



+l3Ej / ^ x^j sin( 




(45) 



where with Xi denoting the x coordinate 

of the ith grain in imits of the lattice constant. The 
universal jump condition of the helicity modulus p7| 



T = 2/'k(3Ej = 2TEo/ttEj 



(46) 



then determines the critical coupling (Ej / Eq)c, the val- 
ues of which are displayed in Table [II, for Ci/Cq — 
and 1 at various temperatures. Table III also displays 
the results from the SCA replacing the Kronecker delta 
by its average, which has yielded Figs. || and |o[ The 
two results indeed exhibit behaviors qualitatively in ac- 
cord with each other. In particular, for Ci/Cq — 0, both 
results give the value of {Ej /Eq)c larger at T 0.1 than 
at r = 0.2, strongly suggesting the reentrant transition. 

From a quantitative viewpoint, on the other hand, the 
agreement between the two results is not satisfactory, 
especially at low temperatures, which seems to indicate 
that the SCA fails at low temperatures. To examine this, 
we have also measured in the MC simulations the aver- 
age value of 5ni,nj for Ci = and display the results 



in Fig. [T|. At high temperatures, the MC data indeed 
agree well with the results from Eq. (41) which is rep- 
resented by the solid line in Fig. |ll|. However, at low 
temperatures the difference is apparent unless Ej/Eq is 
sufficiently small. In particular, for Ej > Eq, the MC 
simulations give the value approaching unity as the tem- 
perature is decreased to zero, whereas the value from 
Eq. (^) approaches zero. This is the case that the sec- 
ond term in Eq. (44) is dominant over the first term, 
which leads to the ordering of both {9i} and {rii}. This 
reflects strong correlations between the two, which have 
not been taken into account in the SCA. It is thus con- 
cluded that the quantitative behaviors of the phase di- 
agrams in Fig. P are accurate only at sufficiently high 
temperatures as already shown in Table 
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V. CONCLUSION 

We have investigated the quantum phase transitions 
in two-dimensional superconducting arrays with various 
values of external charge frustration at both zero and 
finite temperatures. We have considered general capaci- 
tance matrices, allowing both the self- and the junction 
capacitances. In the absence of dissipation, the charge 
variable in the array changes in discrete quanta of 2e, 
and the phase variable becomes compact in the inter- 
val [— tTjTt). We have studied in detail zero-temperature 
phase transitions at various values of the charge density 
and also considered the system at finite temperatures 
with and without the Josephson energy. 

The states of the system at zero temperature have 
been determined from the investigation of the crossing 
of the ground-state level and the first-excited-state level. 
At vanishing charge density ((g^) = 0), we have consid- 
ered the commensurate insulating state, where qi = 
at all sites, and two types of the incommensurate super- 
conducting state, i.e., the single-charge (SC) type and 
the charge-dipole (CD) type. Unless the self-capacitance 
vanishes, the SC type incommensurate state has been 
found to have lower energy than the CD type one except 
at very small values of the charge frustration. From the 
comparison of the energies of the commensurate insulat- 
ing state and the SC type superconducting state, com- 
puted via the second-order perturbation expansions, the 
phase boundary has been obtained. The validity of the 
second-order expansion has been checked for the array 
in the self-charging limit, which reveals that the inclu- 
sion of the third-order term does not change significantly 
the phase diagram. At finite charge densities, identifica- 
tion of the commensurate state and of the first-excited 
incommensurate state is a highly nontrivial problem. We 
have thus adopted the simulated annealing method and 
applied the second-order perturbation theory to obtain a 
number of lobe-like insulating phases such as (qi) — 1/2, 
1/3, 1/4, and 1/5 in 10x10 and 12x12 arrays. As the sys- 
tem size is increased, more insulating lobes are expected 
to be observed. The effects of the junction capacitance 
have been investigated in detail in the insulating phases 
of (qi) = and 1/2: As Ci/Cq is increased, the insulating 
region shrinks monotonically for {qi) ~ 0, while the in- 
sulating region for {qi) = 1/2 expands with Ci and then 
shrinks, with its maximum area at Ci/Cq ~ 0.1. The 
MC simulations have displayed that the lowest excited 
state changes as the junction capacitance is increased, 
which in turn indicates that the existence of the super- 
solid phase with {qi) = 1/2 depends on the interaction 
range of charges: The supersolid phase exists for suffi- 
ciently large values of the ratio of the junction capaci- 
tance to the self-capacitance. In the supersolid phase, it 
has also been shown that the point defect in the charge 
configuration plays an important role. 

At finite temperatures, thermal fiuctuations as well 
as quantum fluctuations play an important role in the 



phase transitions of the system. In the limit of vanish- 
ing Josephson coupling strength, we have performed the 
MC simulations to obtain {qi) as functions of the charge 
frustration at various temperatures, and pointed out the 
possibility of the devil's staircase structure in the thermo- 
dynamic limit. As the temperature is increased, higher- 
order insulating states disappear but some lower-order 
ones still remain at sufficiently low temperatures. In the 
presence of the Josephson coupling, we have obtained the 
effective classical Hamiltonian via a variational method 
in the path integral formalism. A self-consistent approxi- 
mation has been applied to the resulting Hamiltonian and 
phase boundaries have been computed for Ci/Cq = and 
1. A reentrant transition a.i q^ — Q has been observed for 
the array in the self-charging limit (Ci = 0), while an in- 
sulating lobe-like structure has been found near qx = 0.5 
in the case Ci/Cq = 1.0. The validity of the formalism 
used here, the self-consistent approximation (SCA) on 
the effective classical Hamiltonian, has been checked us- 
ing MC simulations in the absence of charge frustration. 
Although the reentrant behavior of the array in the self- 
charging limit has also been observed in MC simulations, 
the SCA appears to be accurate only at sufficiently high 
temperatures. 

The quantum phase transitions of the superconduct- 
ing arrays have been studied theoretically by many au- 
thors via mean-field approximations such as the coarse- 
graining approximation ||^ and via the semiclassical 
methods including the variational methods [p|jic|] and the 
WKB approximation Although the mean-field ap- 
proximations are in general expected to give better re- 
sults at zero temperature, where the dimension of the 
system is effectively increased by one, than at finite tem- 
peratures, it is well known that the coarse-graining ap- 
proximation gives unphysical results even at zero tem- 
perature, for the arrays without self-capacitances. On 
the other hand, the semiclassical methods should be ap- 
plied with great care at zero temperature, where quan- 
tum fluctuations are strong. In comparison with these 
methods, the formalism adopted here to study the zero- 
temperature phase transitions, perturbation expansion 
combined with simulated annealing, appears to be more 
reliable and can be applied to other systems systemati- 
cally. 
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FIG. 1. Schematic diagram of the energy levels En versus 
the control parameter K. In (a) the crossing point Ki of the 
ground state and the first excited state corresponds to the 
transition point. In (b), on the other hand, the true transi- 
tion point is given by K2, at which the ground state meets 
the second excited state. 
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FIG. 2. Boundary for the first excited state when (qi) = 0. 
In the region enclosed by the curve the creation of a charge 
dipole (CD) at one bond gives the lowest excitation, and the 
creation of a single charge (SC) at one site is the lowest ex- 
citation in the region outside the curve. For Ci/Co < 14.116 
the lowest excitation is SC irrespective of Qx, and in the near- 
est-neighbor charging limit (Co = 0) the lowest excitation is 
CD for Qx = and SC for Qx 0, respectively. 
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log(Ci/Co) 

FIG. 3. Zero-temperature phase diagram for {qi) = 0. The 
system is insulating in the region below the surface, which 
shrinks monotonically as Ci/Co is increased. 




FIG. 4. Zero-temperature phase boundaries for (qi) = in 
the self-charging limit, computed via the perturbation expan- 
sion up to 0{Ej"). The dotted line, the dashed line, and 
the solid line correspond to the first-order (n = 1), the sec- 
ond-order (n = 2, and the third-order (n = 3) calculations, 
respectively. 




FIG. 5. (a) The commensurate charge configuration with 
(qi) = 1/2 at zero temperature, where the empty and the filled 
circles denote the vacant {qi = 0) and the occupied (by single 
positive charges; qi = 1) sites, respectively. This commensu- 
rate state is found to be destroyed by a single negative charge 
excitation (q^ = —1) for Ci/Co = 1.0 and 5.0, as shown in 
(b), where the location of the point defect is indicated by a 
small square. When d/Co=0.1, on the other hand, the low- 
est excitation is found to have the configuration (c), which 
has qe = —12. 
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FIG. 7. Zero-temperature phase diagram for (qi) = 1/2, 
where only single negative charge excitation (g^ = — 1) has 
been considered. In the region above the surface, charges 
form a 2 X 2 superlattice, as shown in Fig. ^ (a), leading to 
the insulating phase. As Ci/Co is increased, the insulating 
region expands at first ( for Ci/Co ^ 0.1) and then shrinks. 
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FIG. 6. The zero-temperature phase boundaries for Ci/Co 
= (a) 0.1, (b) 1.0, and (c) 5.0, obtained via the second-order 
perturbation expansion, with the help of the simulated an- 
nealing method to find the charge configuration. Various 
lobe-like structures, inside of which the array is in the in- 
sulating phase, are found in both the 10 x 10 and the 12 x 12 
array. As we increase the array size, additional insulating 
lobes are expected to be observed. As Ci/Co is increased, 
the insulating lobes become narrower, with the central posi- 
tion approaching (qi) = qx- In the inset of (a), the insulating 
lobes with (qi) = 1/3, 1/4, and 1/5 are shown. 
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FIG. 8. Monte Carlo results of ((j'i)T as functions 
of for arrays of sizes L = 10 (left) and 12 (right) 
with = and Ci/Co = (a) 0.1, (b) 1.0, and 

(c) 5.0. Prom the left to the right, the curves corre- 
spond to temperatures T = 0.0,0.006,0.01,0.05, and 0.1 
in (a), T = 0.0,0.01,0.014,0.04, and 0.1 in (b), and 
T = 0.0,0.004,0.008,0.02, and 0.05 in (c), respectively. Here 
the temperature has been measured in units of i?o/fes, and 
all the curves have been shifted in the horizontal direction for 
clarity. More steps are expected to be observed in systems of 
larger sizes. 
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FIG. 9. Phase diagrams in the self-charging limit (Ci = 0) 
at temperatures T = 0.1, 0.2, 0.5, 0.7, and 1.0. Supercon- 
ducting region (the right-hand side of each curve) is shown to 
expand as qx is increased. Reentrant behavior is observed at 
Qx = 0. 
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FIG. 10. Phase diagrams for Ci/Co = 1.0 at temperatures 
T = 0.1, 0.2, 0.3, and 0.4. At T = 0.1, there exists a lobe-like 
structure near qx = 0.5. 
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TABLE I. Values of qe in the lowest excited states of the 
arrays with size (a) L = 10 and (b) L = 12. It is shown that 
\qe\ tends to decrease as Ci is increased. 
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FIG. 11. Monte Carlo results oi 5 = ((5„,,„^.) for 
Ej/Eq = 0.1, 0.6, and 1.0, ai = in the self-charging 
limit (Ci — 0). Results from the self-consistent approxima- 
tion (SCA), which are independent of Ej/Eq, are also shown 
by the solid line. When the temperature is sufficiently high, 
the SCA is shown to give reliable results; at low tempera- 
tures, SCA results are in agreement with the MC results only 
for sufficiently small Ej/Eq. Dotted lines are merely guides 
to the eye. 
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TABLE II. Step width for (qi) — 0, 1/3, and 1/2, in the L = 12 array at zero temperature, obtained from the data in Fig. 
The numbers in parentheses denote the maximum error in the last digits. 



(90 


Ci/Co = 0.1 


Ci/Co = 1.0 


Ci/Co = 5.0 





< < 0.3650(5) 


<qa: < 0.1275(5) 


< g^ < 0.0395(5) 


1/3 


0.3900(5) <qa: < 0.4055(5) 


0.314(1) < < 0.3800(5) 


0.317(1) < g, < 0.356(1) 


1/2 


0.4130(5) < < 1/2 


0.4290(5) < g^ < 1/2 


0.4725(5) < go, < 1/2 



TABLE III. The critical values of Ej/Eq at g^ = for C\/Cq =0 and 1.0. The results from Monte Carlo simulations and 
those from Figs. ^ and |l^ at temperatures T — 0.1, 0.2, and 0.5 are displayed in the left and the right columns, respectively. 
Both results display reentrance for C\/Cq — 0. The numbers in parentheses denote the maximum error in the last digits. 



T 


(Ej/Eq), for Ci/Co = 


(Ej/Eo)c forCi/Co = l 


0.1 


0.642(1) 


0.88(1) 


0.158(2) 


0.28(1) 


0.2 


0.558(2) 


0.81(1) 


0.242(4) 


0.34(1) 


0.5 


0.766(2) 


0.91(1) 


0.53(1) 


0.57(1) 
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